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Our progress over the last year has been along several dimensions. 

a. Exploration and understanding of Earth Observatory System (EOS) mission 
with available data from NASA 

b. Comprehensive review of state of the art techniques and uncovering of 
limitations to be investigated (e.g. computational, algorithmic ...) 

c. Preliminary development of resolution enhancement algorithms 

1) Background 

• Radiometry 

With the advent of well-collaborated satellite microwave radiometers, it is now possible to 
obtain long time series of geophysical parameters that are important for studying the global 
hydrologic cycle and earth radiation budget. 

Over the world’s ocean, these radiometers simultaneously measure profiles of air temperature 
and the three phases of atmospheric water (vapor, liquid, and ice). In addition, surface 
parameters such as the near surface wind speed, the sea surface temperature, and the sea ice 
type and concentration can be retrieved. 

The special sensor microwaves imager SSM/I has wide application in atmospheric remote 
sensing over the ocean and provide essential inputs to numerical weather-prediction models. 
SSM/I data has also been used for land and ice studies, including snow cover classification 
measurements of soil and plant moisture contents, atmospheric moisture over land, land 
surface temperature and mapping polar ice. 

The brightness temperature observed by SSM/I is function of the effective brightness 
temperature of the earth’s surface and the emission scattering and attenuation of the 
atmosphere. Advanced Microwave Scanning Radiometer (AMSR) is a new instrument that 
will measure the earth radiation over the spectral range from 7 to 90 GHz. Over the world’s 
ocean, it will be possible to retrieve the four important geographical parameters SST, wind 
speed, vertically integrated water vapor, vertically integrated cloud liquid water L. 


• Resolution Enhancement 

Due to the diffraction of differences in antenna pointing directions, microwave radiometer 
antenna gain patterns of different frequency may differ in size or location on the earth surface 
while physical inversion algorithms may assume that observations describe consistent 
locations. Therefore, the mismatch in the resolution becomes a critical problem when 
observations at many different frequencies must be combined to retrieve single parameters, 
which requires that the area on the surface imaged by different channels be the same. 

The obvious solution to this problem is to average the high-resolution measurements to 
match the lower resolution data. However, this solution is not recommended due to the non- 
linearity in the physical model that relates the measured brightness temperature to the 
geophysical parameters. A much more desirable solution, is to enhance the low- resolution 
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data to match the higher resolution data. The theory is based on the fact that the density of 
the measurements made by the space borne radiometer is higher than the resolution of the 
instrument, which means that it is possible to take advantage of over sampling to reconstruct 
a higher resolution image from low-resolution data. 


2) Survey of Image Restoration Techniques 

The purpose of image restoration is usually formulated as the estimation of an improved 
image 5 of the original image 5 when a noisy blurred version Z given by: 

Z=HS+V (1) 

is observed, where the blurring matrix H is known and some statistical knowledge about V 
and S is assumed to be available. Specifically, we will assume that the image model 
S=AS+W is feasible. The direct inversion of the matrix H does not lead to useful restoration 
because of the ill conditionedness of H. 

Procedures to stabilize the inversion of an ill conditioned matrix are called regularization 
methods and make nearly always use of a priori knowledge about the original image and the 
noise. 

The main objective in solving ill posed problems is the construction of physically acceptable 
and meaningful approximation of the true solution that is sufficiently stable from a 
computational point of view. If we are to obtain a useful approximate solution to (1), we 
must modify the problem in such a way that the solution becomes less sensitive to noise in 
the observed image (well posed problem). At the same time the solution of the modified 
problem must be close to the solution of the original problem. The process of achieving a 
compromise between these conflicting goals is referred to as regularization and is typically 
controlled by one or more regularization parameter. Nearly, all concepts used in 
regularization are based on incorporating a priori knowledge about either the true solution or 
the noise into the algorithm which solves the image restoration problem 


2-1) Tikhonov-Miller Regularization: 

The idea is to define a criterion to select an approximate solution from a set of feasible 
solutions. On the basis of the observed model (1) it is plausible that the class of feasible 
solutions is described by : 

0(5) = 1 Z - HS\\ < |F| = e (2) 

The bound e is related to uncertainty or noise in the observed image Z and can usually be 
estimated from smooth image region. Tikhonov defined the regularized solution as the' one 
which minimizes the stabilizing functional 0(5) in the set of feasible solutions. Although a 
wide class of different stabilizing functionals is available, including for example the max 
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power and max entropy measures, usually a stabilizing functional of the following form 
chosen to facilitate the mathematical analysis of the problem 

Q(S) = \\CS\\ (3) 

Here C is a real valued matrix of size MxN known as the regularizing operator. The 
computation of the regularized solution now reduces to the minimization of (2) subject to (3). 
Using the method of undetermined Lagrange multiplier we need to minimize the objective 
function: 


0(5) = \Z - //S| 2 + a||C5|| 2 (4) 

with respect to S. The regularization parameter «is chosen so that (2) is satisfied with 
equality. [1] proposed a empirical criterion based on CRESO (composite residual and 
smoothing operator) for finding the optimum value of a"'" which gives the solution S closest 
to its initial value. According to this criterion a^'is the value a that gives the first local 
maximum of the function: 

C(*H|S| ! + 2«j;M 2 (5) 


Among the solutions satisfying a reasonable choice is the one which minimizes 0(5) called 
the Tikhonov Miller regularized solution. This minimization is straight forward and leads to 
the following solution: 

(H , H + aC'Cy i S = H‘Z (6) 

provided that (H'H + aC'Cy' is invertible. 

Iterative Solutions 

Iterative solutions offer the advantage that no matrix inverses need to be implemented and 
that additional deterministic constraints can be incorporated into the solution algorithms. 

To derive the iterative Tikhonov Miller regularized restoration algorithm, many methods 
were proposed [2]. 

• Steepest Descent Method: 

The steepest descent method is used to minimize the objective function 0(5) by using the 
following iteration: 

5* +l = 5* + (3r k 

= S k - (3({H'H + aC‘C)S k - HZ) 
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= (/ - a/3C'C)S k + PH' (Z - HS k ) 


(7) 


This iteration reduces the to Van Cittert iteration if a =0 which is the case of no 
regularization. 


S k " =S k +0H‘ (Z- HS k ) (8) 

The term (I-aflC’C) in () behaves like a low pass filter suppressing the noise in the 
iterates. The characteristics of this term are obviously related to the properties of the original 
image because the regularizing operator C is closely related to the image model A. 


Analysis of Convergence Speed 

The convergence speed of an iterative algorithm towards its limiting solution can be 
quantified by its convergence rate. An iterative scheme converges geometrically with order R 

if the error | S'* - S’"! for k sufficiently large can be given as 

||S* -5*||= K\\S k -S c0 ||\/?>1 (9) 

A large value of K corresponds to a slow convergence while a small value indicates fast 
convergence. Further, the larger the convergence order is, the faster the algorithm converges. 
For R=T , the process is said to converge linearly in which case (9) is conveniently written as: 

||S* -S"||< K k \S° -S' 0 !* (10) 

and 0<££1. For K=0, the convergence is said to be superlinear or in other words the 
iterative process terminates within a finite number of iterations. While for K=l, the 
convergence is sublinear. 

Iterative schemes which minimizes a quadratic objective function O(S) by the method of 
steepest descent can all be regarded as members of the following generic iterations 

S k+] =S k +/3((h~BS k ) (11) 

for nonsingular matrix B, we have S K = B 'H . The following relation can then be 
established: 

S* +l -S” = (I- pB)S k +J3h-S 00 

= (I-/3B)S k + pBS x - 5" 

= (I-/JB)(S k -S") (12) 
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I S k -S"|<||/-,®||||s /t -S*f (13) 

Here || I - j3B\ denotes the norm of the matrix which is given by its largest eigen value. Thus, 
if we let p mm be the eigen value of B with smallest norm then (13) becomes : 

I?* - 5 ,cc | < |l - (3p mm |||<S* - S"|| ff (14) 

Hence, iterative schemes of the type (12) converges linearly with a convergence rate of 
|1 - pp mm \ . For the iteration (12) , we have (/ - (3B) = / - fi(H' H + aC'C ) and p mm is the 

smallest eigen value of (H' H + aC'C ) . 

Since the smallest eigen value is nearly always close to zero in an image restoration 
application, the iteration (12) converges slowly. 

In general, the convergence rate depends on the regularization parameter a because p mm is a 
function of a . If a is increased, i.e. if the restoration problem is regularized more strongly, 
p mm attains a larger value by which the speed of convergence increases as well. Also, /? was 
assumed to have fixed value satisfying the condition: 

0</?< (15) 

J/Tnax I 

where p max is the largest eigenvalue of the matrix ( H'H + aC'C) . 

But since /? controls the convergence rate as follows from 0, it is desirable to optimize its 
value at each iteration step. 


• Conjugate Gradient Method 

Motivated by the desire to achieve more rapid convergence, the method of conjugate 
gradient has been successfully used in optimization theory. Conjugate direction methods 
which were originally introduced for purely quadratic problems, can be viewed as special 
orthogonal expansion of the solution of the minimization problem. This expansion is 
generated by making use of information from previous iteration steps. One advantage of this 
method is its convergence in a finite number of iterations steps when exact arithmetic is used 
(no rounding errors) the convergence is superlinear. When nonexact arithmetic is used or the 
problem is non quadratic, the method will no longer converge in a finite number of iterations 
because the conjugacy condition will no longer holds. 

It has been experimentally shown, however, that the conjugate gradient method converges 
always faster than the method of steepest descent, while the computational complexity is not 
significantly increased. The basic form of the conjugate gradient algorithm which thus 
represents an alternative to the iteration (12) is given by: 
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(16) 


r* = - - V S <D(5) | s = -(//'// + aC'C)S k + H'Z 
2 * 


Pk = r k +YkPk 
S M =S k +f3 kPk 


In this scheme, S k is modified in the direction of the vector p k instead of the steepest 
descent direction r k . The parameter y k controls the conjugacy of the subsequent direction 
p k . For y k — ¥ 0 , the iterations reduce to (7). It can be shown that the values y k and are 
given by: 


Yk = 



(17) 


P k 


rjPk 

\\H Pk ( +a\\Cp k ( 


(18) 


The conjugate gradient method converges linearly with a convergence rate given by: 

r . , x* 

| S k -5®||< K 

Although (19) is rather upper bound to the convergence rate of the conjugate gradients 
iterations (for example, it does not show that the method of conjugate gradient converges 
within a finite number of iterations, it is still has a higher convergence speed than the method 
of steepest descent. In addition to this, experimental results exhibit a convergence speed 
much larger than indicated by the bound (15). 

Observe that the difference in computational complexity between the method of steepest 
descent and the conjugate gradients algorithm merely consists of the computations of p k and 
y k which is insignificant compared with the other computations required within a single 
iteration step. 


\1 Anax | Y I Ai 

max \l\ P m 


(19) 


• Preconditioned Conjugate Gradient Methods 

The convergence rate of iterative methods depends on the spectrum of matrix H. Generally 
speaking, the method converges faster if H has smaller condition number or clustered eigen 
values. In order to accelerate the convergence rate we may consider the solution of a new 
equivalent system of equations: 


P-'Z = P~ ] HS (20) 
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instead of the original one where P can be easily interpreted and P~ X H has better spectral 
property than H. the matrix P is usually called the preconditioner of the matrix H. To 
impalement preconditioning, the matrix vector product P~'v for a certain vector v is 
performed at each iteration. Various preconditioning techniques have been introduced in 
solving elliptic partial differential equations and Toeplitz systems of equations. 

Preconditioned CGN method was introduced in [4], where the conjugate gradient methods is 
applied to the normal equations H'Z = H'HS when the matrix H is not symmetric positive 
definite since H' H is symmetric. 

Consider the model: 


Y=AS (21) 

Where A is a block Toeplitz matrix (symmetric positive semi definite). A good 
preconditioner for A is a matrix P that approximates A well in the sense that the spectrum of 
the preconditioned matrix P~'A 

is clustered around 1 or has small condition number and that can be inverted efficiently. 

One preconditioning technique has been generalized from point to block Toeplitz matrices by 
Ku and Kuo [3], They proposed to use a block circulant matrix as a preconditioner for the 
block Toeplitz matrix since the block circulant matrix can be easily inverted via 2-D FFT. 

Let A be a block Toeplitz matrix containing NxN blocks with MxM elements per block: 


A 


A= 


A 




* 2 -N 


*1 -N 


(22) 


^V-l ^V-2 


A 


Where A„ with |«| < V - 1 are MxM Toeplitz matrices with elements [A n ] = a : _ jn , 
1 <i,j<M 

We use the following MNxMN block circulant matrix as its preconditioner 


Pn-i 


Pi 


P = 


N~\ 


*jV-2 


Where 

r 

p 


p + p 


N-n* 


n = 0 

1 <n < N -n 


(23) 
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and where P„ is circulant matrix with a 0/l ,a_^ +a u _ w a_ 2/l +a w _ 2 „ +a x _ Mn +a^ as the 

first row. 

From the construction, we see that P is an approximation of A and the matrix vector product 
P~'v for any v can be computed via 2-D FFT. 

The PCGN method is summarized as follows: 


Arbitrary x 0 , r 0 = y - Ax 0 , p 0 = 0, /? 0 = 0 


Iteration k=l,2,.... 


z k-\ ~ A'r k _\ 

P k = z k - 1 + PkPk-\ 
= (z» ± A d ) 

*" Hftf 

X k — X k~ 1 + &kPk 

r k = r k-\~ a k A Pk 

fa. A) 

{z k _„A'r k _ x ) 


(24) 


where M is the preconditioner. M = P'P . 

The convergence rate of the preconditioned iterative method depends on the eigen value 
distribution of the preconditioned matrix (P'P)' 1 . 

The preconditioned iterative methods converge faster if ( P'Py'A'A has clustered eigen 
values and /or small condition number. 

Kuo [3] studied the eigenvalues distribution of the matrices A' A and ( P'Py'A'A with 
different values of the image size where they showed that the eigen values of the 
preconditioned matrices have much better clustering property that the original ones. 


2-2) Algebraic reconstruction techniques ART 

ART is a reconstruction technique that was recognized to be identical to the algorithms 
of solving the system of liner equations Z=HS where Z,H,S are as defined in (1). 

ART are based on an iterative process, which starts from an initial approximation SO to the 
image vector. In an iterative step, the current iterate/?* is refined (or corrected) to a new 
iterate S k+] by taking into account only a single measurement, say the ith measurement and 
changing only the image values of the pixels which intersect this measurement. The 

discrepancy between the measurement Z, and the pseudo projection data Ys h v S j obtained 

1 
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from the current images'* is redistributed among the pixels along the ith measurement 
proportionally to their weights A, 7 in the whole measurement. 

In this way, the pixel values covered by the ith measurement are corrected to conform with 
the ith measurement without changing the rest of the image. Denoting h 1 = (A y )" =l as a vector 

in R n . The following algorithm describes this process: 

S° e R m 
Typical step 

S* +1 = S k + Z| ~ <hl ;‘ S * - h' (25) 

|h‘|| 

X " 1 ■ II II 2 

Where: < w,v >= an d ||«|| =<u,u>. 

7=1 

and the measurements are chosen cyclically: i - i k = £(mod M ) + 1 . 

Censor [5] investigate the convergence behavior of ART when applied to ill posed and 
inconsistent systems of equations and proposed to use relaxation parameters which are a 
sequence (^)” =0 of real numbers usually confined to an interval 

£■, < X k < 2 - e 1 , £, , s 1 > 0 


Which appear in the typical step of ART as: 


S 


*+i 




Z,- < h',S* > 



h' 


(26) 


The relaxation parameters allow are to overdo or undo the orthogonal projection prescribed 
by ART. 


2-3) Block ART 

Block ART were first introduced by Eggermont et-al [] to solve the system of ill posed linear 
equations Z=HS, where the MxN matrix H may be partitioned in two ways as: 


'K~ 


i 

i 

K 

— 


_K*_ 


"" i 

i 


here h, is an N dimensional vector and its transpose h* constitutes the ith row of H. each 
matrix H, is an LxN block of rows of H and assuming there are K blocks of width L 
obviously means that M=LK. 
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The M dimensional vector Z is partitioned as 

Z' =(Z„Z 2 , ,z lk ) = l z;,z',...,z' K ) 

where z, is a block of length L of elements of Z. 


The block ART is the following group of iterative variations 

Initialization 

S° e R N is arbitrary 

Typical step: 

S k+ ' = S k + H‘Q k (Z, - H,S k ) (27) 


where Q* is an LxL matrix called relaxation matrix, H’ is the Moore-Penrose inverse of 
H, and the blocks with respect to which the iteration is performed are chosen in a cyclic 
manner: / = i k = &(mod M) + 1 , k= 1,2,... 

The case Q* = Z k I with I standing for an LxL identity matrix coincides with ART 
algorithm when L=l. 

In addition, Long et-al, showed that block AART is equivalent to a least squares estimate 
in the limit of infinite iterations based on the minimization problem: 


Minimize ||S|| 2 subject to Z=HS. 

Which is equivalent to Tikhonov regularized solution. 


2-4) Entropy Regularization (MART) 

The use of entropy is rigorously founded in several areas. In image reconstruction 
from projections several authors advocate the maximum entropy approach. From the 
standpoint of information theory, this approach is conceptually attractive. It yields the 
image with the lowest information content consistent with the available data. 

Thus with this approach one avoids introducing extraneous information or artificial 
structure. The problem of reconstructing a source from a finite number of views is known 
to be indeterminate. A maximum entropy method thus seems attractive for this problem 
especially when the available projection data are incomplete or degraded by noise errors. 
By Frieden’s model [6], the most likely object scene implied by given image data is 
found to obey a principle of maximum entropy. 

Entropy optimization refers to the problem of maximizing the functional 

n 

f (S) = S / In over the constraint Z=HS. 

./=! 
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( 28 ) 


A typical problem would be 

n 

Maximize: - lnS, subject to Z=HS 

j = i 

and S' > 0 

A variety of iterative algorithms are available to solve () but no practical direct form 
solutions are known. 

Lent [5] suggested the following algorithm and proved that it actually converges under 
some reasonable conditions to the solution of the maximum entropy problem. 


= S (- 


-y 


<h',S k > 
where X k is a relaxation parameter. 
Equations are taken cyclically, i.e. i = i k 


(28) 


= £(mod M ) + 1 . 


2-5) Variational Methods 

In the previous sections we showed that the classic setting for solving the image 
reconstruction inverse problem is the least squares minimization. However, regularization 
is needed when solving the ill posed problems to neutralize some of the ill conditioning 
but it also removes sharp edges and similar distinguishable features. Hence, the choice of 
the regularizing function is an issue. 

In variational approach, an image is defined as a real function S: S' : £2 — > R , where Q is 
an open bounded subset of R 2 . 

The formulation of the image reconstruction using variational methods is based on 
solving the optimization problem 
Minimize O(S') subject to Z=HS 
i.e. 

0(S) + y||Z-//S|| 2 j (29) 

The critical issue is the choice of the variational integral. One of the widely used 
variational integrals based on geometrical argument is the total variation proposed by 
Ruden et-al and defined as follows 

7V(S) = J|VS|c/<y (30) 

Q 

the total variation method basically consists of finding an estimate S for the original 
image S with the smallest total variation among all the images satisfying the noise 
constraint \\Z-HS\\ = cr 2 where a 2 is assumed known. Roughly TV permits image 

intensity to have sharp jumps, but it limits spurious oscillations. Regularized least 
squares, on the other hand, tends to smooth out sharp jumps because it controls the 
second derivative of the image intensity. 


S - arg minj 
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An iterative solution for the optimization problem can be obtained using gradient descent 
method as follows: 

S* +l =S k +a[v(VS*)-A(Z-//S*)] (31) 

Ben Hamza and Krim [7] introduced an entropic variational approach based on the 
concept of negentropy variational integral which is defined as: 


H(S) = J H(\VS\)dco = J|VSj log|' VS\dco (32) 

n a 

Then the minimization problem can be written as: 


S = arg mini J|VSj log(|VS| + — ||Z - HSf dco , 

s U ^ J 

and the iterative solution can be written as: 


-.i+l 




1 + log VS* t 

V( 3-t — VS* ) - \(Z - HS k ) 

VS* 


(33) 


(34) 


3) Spatial Resolution Enhancement of SSM/I Data 

The special sensor microwave imager SSM/I, a satellite passive microwave 
radiometer observes a large portion of the earth’s surface and measures the top 
atmosphere brightness temperature at 19.35,37,85 GHz in both the horizontal and the 
vertical polarizations and at 22.235 GHz channel in the vertical polarization only. 

The instrument flies at an altitude of 833km and scans the earth in a conical scan 1400km 
wide. The SSM/I samples data in seven channels covering four microwave frequencies - 
table (1). The lower frequencies of the instrument (19,22,37) are sampled every 25km 
along the scan and every 25km along the track, while the other high frequency is sampled 
12.5km. 

Due to the physical limitations on the SSM/I antenna size the data have relatively poor 
spatial resolution. The spatial resolution of a passive microwave sensor depends both on 
the antenna size and frequency. Because the SSM/I instrument uses on parabolic reflector 
for all four frequencies, the spatial resolution improves with increasing frequency. 

To retrieve geophysical parameters from SSM/I measurements take at different 
frequencies requires that the area on the surface imaged by different channels be the same 
in order to obtain more accurate measurements of geophysical parameters, it would be 
preferable to improve the resolution of data at low frequencies to that of higher 
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frequencies. Table(l) shows the SSM/I channels and the corresponding sizes of foot 
prints and sampling intervals. Correction techniques for matching the resolution of all 
frequencies of the SSM/I to the 25 km spatial resolution of the 37 GHz were presented in 
[8]-[ 13]. 


Table 1 


Channel 

Frequency 

Polarization 

Resolution (km) along 
Track Scan 

Sampling 
interval (km) 
Along 
track/scan 

19.35 

V/H 

69 

43 

25 

22.235 

V 

60 

40 

25 

37 

V/H 

37 

28 

25 

85.5 

V/H 

15 

13 

12.5 


Backus & Gilbert first proposed the methodology for the correction techniques. In 
essence, the theory contends that if the density of the measurements made by the satellite 
is higher than the resolution of the instrument, then it is possible to find a linear 
combination of surrounding measurements that yields results at a higher spatial resolution 
than the original data. The tradeoff of r higher spatial resolution in this procedure is the 
rapid amplification of instrument noise. Other techniques based on algebraic 
reconstruction techniques provide an estimate of the brightness temperature of each 
element of a rectilinear grid of pixels. The problem reduces to the solution of a system of 
linear equations that is ill conditioned by using iterative algorithms. ART requires less 
computations than BGI methods, but on the other hand, BGI methods provide a tradeoff 
between noise and resolution enhancement. 


3-1) BGI Algorithm 

An SSM/I measurement can be modeled as a product of the surface brightness 
temperature and the antenna pattern. The ith measurement T A (i) is obtained by 
integrating the product of the surface brightness T B (x ,y ) and the antenna at the surface 

G,(x,y) 


T A (i) = G~ 1 JJ G,(x,y)T b (x,y)dxdy (35) 

G,-' = JJ G,(x,y)dxdy (36) 
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Where integrals over surface area are corresponding to he non-negligible gain of the 
antenna. The dependence of G on I arises from bore sight pointing of the antenna, which 
changes as the antenna scans the surface. 

The antenna pattern acts as a low pass filter of the surface brightness limiting the 
effective resolution of the measurements. The BGI algorithm is an inversion method for 
solving the integrals equations. The algorithm is used to determine surface brightness 
from integrated over lapping antenna patterns. When employed for spatial resolution 
enhancement, the BGI method produces a weighted least squares estimate of the surface 
brightness temperature on a rectilinear surface grid finer than the intrinsic resolution of 
the sensor. Given a set of antenna temperature measurements T A (i) with associated 
antenna gain patterns G, (x, y) , the algorithm estimates the brightness temperature 
T B (x , y ) for each pixel (jc j , y J ) of the enhanced resolution image. 

To estimate Tb for a given pixel, the BGI method uses linear combination of N nearby 
measurements 

M 

< 37 > 

coefficients a tJ are determined from measurement geometry and the noise correlation 

matrix. These coefficients are different for every pixel due to the varying antenna 
geometry over the swath. 



Figure 1: Reconstructing high-resolution brightness temperature from T A (i) measurements. 


Substituting (1) into (3) 


f \GX*,y) T ,(x,y) = Jjl [a,G(x,y)\,(x,y)dxdy 

1= 1 


(38) 
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If one can find coefficients a tJ such that the bracketed term becomes the dirac delta 
function 5{Xj - x,y J - y) , then (38) would be an exact solution to the problem and (37) 

could be used to find the actual scene brightness temperature at any desired resolution. In 
practice, however, the finite number of measurements makes it impossible to find such a 
set of coefficients a v . 

Instead, we confine ourselves to selecting a set of coefficients a i} which causes the 
bracketed term to most closely produce some desired behavior. In particular, we look for 
a set of coefficients a which make the bracketed term appear like the gain function 

G(x,y ) of the 37 GHz channel. It must be understood that T B (Xj,yj ) is no longer the 

true scene brightness temperature but the brightness temperature as it would be seen at 
the 37 GHz resolution. By using the actual 37 GHz antenna function as the match 
function, no correlation is necessary for the 37 GHz s. The 37 GHz antenna function also 
has a footprint which is most similar in size to the 25km spatial sampling distance if the 
instrument. 


To find the coefficients a , consider the quantity: 


l M 

Qo = || Y,a IJ G l (x,y)-G 3 \x,y) 
/-! 


dxdy 


(39) 


Where G 37 ( x,y ) is the gain function of the 37 GHz channel. 

Minimizing Q 0 with respect to a tj subject to the constraint that proper normalization be 
presented. 


r r M 

\&afiXx,y)tkdy = 1 (40) 

j = l 


should produce the best correlation coefficient possible in the least square sense. 

The above treatment assumes the antenna measurements to be exact. In practice, 
however, the antenna temperature T A (i) will not be known exactly, but rather will 

contain some random noise whose variance is given by (A T rms ) 2 

Thus noisy T A (i) can be expressed as 

f A (i) = T A (i) + £, (41) 

Where: T A (i) is the mean value of T A (i) 

£■, : random noise component of zero mean. 

The high-resolution brightness temperature can be expressed as: 

7^ =a'f< =a'T A +a'e (42) 
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where T B is a simplified notation of T^x^y^, a' =[a ]J a 2j ...a Mj ] 
T a =[T a (\),T a (2),...,T a (M)] t . 


And the noise variance of 7^ is: 

a 2 =E[(f Bi -T g ) 2 } 

= E[( ae)(ae)'] 

- £[(aee‘a)] 

= a£[es']a 

= a3a (43) 


Where 3 is the noise covariance matrix of instrument noise process appearing in the 
sampled antenna temperature and is a function of system noise temperature, predetection 
bandwidth and receiver low pass filter when the noise in each sample is uncorrelated and 
equal, then e 2 is proportional to the sum of squares of a tJ times the variance of 
instrument noise. That is: 


(A TJ? 


3 = 


(A7) mi ) 2 J 


(44) 


a diagonal matrix. 

Because we not only want to increase the resolution but rather to avoid introducing any 
drastic increases in the instrument noise. We thus seek a compromise between resolution 
enhancement and noise. 

Rather than Q 0 we minimize: 

Q = Q g cosy + e 2 wsiny (45) 

Where w is a scale factor used to make the two terms on the right hand side of (45) 
dimensionally and numerically compatible, y can be varied between 0 and £_ to place 

various degrees of emphasis on either the resolution enhancement or the noise 
suppression in the estimate of T fl (x 7 ,y ; ) . 

The coefficients a , can now be expressed as: 

a = S' [v cosy- Au] (46) 
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where: 


, _ -1 + uS'Vcosy 

A — ; 

uS u 

S, lk = cos y JJ G,G k dxdy + S IJ (AT rms ) 2 wsin y 
u = [u ] ,u 2 ,...,u M ] r , 

v = [v„v 2 ,...,v w ] r 

u , = JjG^afy, l,k=l,2,...,M 

v, = JJ G 21 G t dxdy 

5 tJ is kronecker delta function, 

Note that these coefficients are different for every pixel due to the varying antenna 
geometry over the swath. 

There are two tuning parameters, the dimensional parameter w and noise parameter y. 
w is arbitrary. Following [1 1] w is set to 0.001. The noise tuning parameter y which can 

TC 

vary from 0 to — controls the tradeoff between the resolution and noise parameter. 

2 

y can be subjectively selected to optimize the resulting image. The optimum value for y 
is dependent on the value of A T used for the noise level. 

In general y is different for each SSM/I channel. [12] developed an objective function 
technique for selecting y for the 19, 22 and 37 GHz channels based on maximizing the 
correlation between the 85 GHz channels and the particular channel of interest. 

Note that when multiple images with exactly the same measurements and pixel locations 
are processed, the BGI enhancement coefficients may be stored and reused. 

The BGI produced image is affected by the definition of nearby N measurements and the 
relative location and gain patterns of the measurements included in the sum (37). 
Restricting the size of the local region defining nearby measurements reduces the 
computational load at the expense of accuracy. Increasing the size of the local are (and 
M) to include additional measurements can improve the accuracy of the resolution 
enhancement but may significantly increase the computational load. 

Previous investigators have used a fixed range of angles around the antenna pointing 
direction to define nearby measurements. This approach can lead to numerical problems 
in the matrix inversion step of the algorithm if the angle range is large enough to include 
directions with very small region. 

Long [14], defined nearby measurements as those measurements that have nonnegligible 
gain at the pixel of interest. A threshold is used to determine if the gain is nonnegligible. 
The measurement is used in (37) only if the relative antenna gain at the pixel of interest is 
greater than the threshold. Setting a lower gain threshold results in more measurements 
being used in () but increases the noise level of the images. 
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The exact number of nearby measurements used at a given pixel is a function of the gain 
pattern and sampling geometry and varies across the swath and for each channel. 

The BGI approach produces adequate enhancement of the spatial resolution to make such 
a correction worthwhile. 

Results obtained in [14] suggests that the BGI technique decreases brightness 
temperature differences by over 50% for the low resolution channels with only a modest 
increase in random noise. 

The correction can also help to better resolve small features which would otherwise be 
lost due to the lack of resolution. This can be important in making better retrievals of 
atmospheric quantities. 

This correction can also be used to average the 85 GHz channels to match the 37 GHz 
resolution when uniform spatial resolution is required. 


(3-2) Enhancement of SSM/I Images Using Image Reconstruction Techniques: 
Problem Formulation and SIR Algorithm 

In order to develop a method for estimating the brightness temperature on a higher 
resolution grid T B (x,y) from antenna temperature measurements T A (/') . 

We consider the desired resolution grid with pixel values T B (x,y ) to be estimated and a 
set of measurements T A (i) whose bounding rectangles are completely contained within 
the region of interest (figure 2). 

Assuming that the brightness temperature is constant within each pixel: 

t a (o = i X Z y^ G > ( 47 ) 

G, x-ly-B 

where 

G, = 'LY.GXx,y) 

and L,R,B,T define a bounding rectangle for the ith measurement. 
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Figure 2: bounding rectangles that approximate antenna temperature measurements 


Denoting the desired resolution brightness temperature value T B (x,y) by vector S where 
S is an N dimensional vector of T g (x,y) in lexicographic order, the radiometer 
measurements T A (i) by vector Z of dimension Mxl and let H be an NxM matrix 
containing the weighting function of corresponding antenna gain function h tJ 
Where 0 <h l] < 1 , which are functions of the antenna gain function and location. 

Then we obtain a matrix equation relating the measurements in Z to H and S as: 

Z=HS+V (48) 

Where V is an M dimensional vector of the noise terms v, , /=/, ...M. 

Depending on the resolution element size and the number of measurements Z , H can be 
either over or under determined. 

The elements of matrix H are computed by determining the intersection of the 
measurements and the resolution elements. We can assume that the entries of H are either 
1 or 0. The resulting matrix is very sparse though it may be extremely large. 

With this formulation, the problem of obtaining estimates of T B is reduced to the solution 
of an ill posed linear system of equations which is similar to classic image reconstruction 
problems in signal and image processing. 


5) Scope of Current and Future Work 

The AMSR instrument was launched in the year 2000 aboard the Japanese ADEOS-II 
platform. AMSR observations will be collected in scans of 196 observations over a width 
of 1600km. The spacing between observations is approximately 10 km for all but the 
highest frequency and approximately 5 km for the highest frequency. Due to the conical 
scan geometry, the density of observations increases at the edges of the scan. Table 2 
summarizes the range of time integrated AMSR footprint sizes. 

The only previous work that has been done to match the resolution of the AMSR 
measurements was that of Ashcroft and Wentz [16], where they used the BGI method to 
match the resolution of AMSR channels. Their work didn’t not include any enhanced 
resolution in order to avoid noise amplification. Instead, each set of actual observations is 
resampled to one of several lower resolutions corresponding to that of 6.9 GHz, 10.7, 
18.7, 36.5 GHz and 89 GHz channels. 

Ashcroft and Wentz [16], have also addressed the artifacts of the BGI method when used 
to match the resolutions of AMSR channels. They analyzed the tradeoff between 
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resolution matching and noise amplification for one combination of AMSR channels at 
one particular point in the scan where they show that because the scan geometry changes 
across the scan, the precise relationship between smoothing, noise and resolution 
matching will also vary. This problem appears in choosing the parameter y which trades 
off resolution and noise as described in the previous sections. 

In addition, they addressed the computationally intensive task of calculating the 
weighting coefficients which is considered as a more challenging problem over the 
SSM/I case because of the higher density of measurements of the AMSR instruments 
which means increased number of calculations. 


Table 2 


Channel 

Frequency 

Resolution (km) along 
Track Scan 

Sampling 
interval (km) 
Along 
track/scan 

6.9 

70 

40 

10 

10.7 

45 

26 

10 

18.6 

27 

17 

10 

23.8 

20 

13 

10 

19.35 

13.5 

43 

10 

36.5 

60 

11 

10 

50.3 

9 

10.5 

10 

52.8 

9 

10.5 

10 

89 

5.3 

5.5 

5 


Motivated by the previous work that has been done for enhancing and matching the 
resolution of SSM/I measurements, we are going to introduce a new algorithm for 
enhancing the resolution of AMSR images by utilizing the techniques described in the 
previous sections taking into account the different nature of the AMSR instrument where 
the number of frequency channels used is increased over SSM/I channels and the 
different sizes of the footprints in addition to the spacing between measurements and 
noise margins introduced by the instrument. 


Our new approach uses to advantage the different resolutions from the antenna by 
constructing a hierarchy of scales which is in turn used to efficiently compute a model 
relating them. This hierarchical structure as shown in Figure 3, lends itself to not only a 
parallelization and simple hardware implementation. When resolutions are harmonized 
across modalities, in a systematic, and reliable way our next objective will be to develop 
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algorithms 


to 


fuse 


data. 
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